Persistence of a Continuous Stochastic Process with Discrete-Time Sampling 
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We introduce the concept of 'discrete-time persistence', which deals with zero-crossings of a con- 
tinuous stochastic process, X(T), measured at discrete times, T = nAT. For a Gaussian Markov 
process with relaxation rate /«, we show that the persistence (no crossing) probability decays as 
[p(a)]' 1 for large n, where a = exp(— fiAT), and we compute p(a) to high precision. We also define 
the concept of 'alternating persistence', which corresponds to a < 0. For a > 1, corresponding to 
motion in an unstable potential (u < 0), there is a nonzero probability of having no zero-crossings 
in infinite time, and we show how to calculate it. 
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Persistence of a continuous stochastic process has gen- 
erated much recent interest in a wide variety of nonequi- 
librium systems including various models of phase order- 
ing kinetics, diffusion, fluctuating interfaces and reaction- 
diffusion processes . Persistence has also been recently 
used in fields as diverse as ecology || and seismology [|| . 
Persistence is simply the probability P(t) that a stochas- 
tic process x(t) does not change sign up to time t. In 
most of the systems mentioned above, P(t) ~ t~ 6 for 
large t, where the persistence exponent 9 is nontrivial. 
Apart from various analytical and numerical results, this 
exponent has also been measured experimentally in sys- 
tems such as breath figures §, liquid crystals §, soap 
bubbles || , and more recently in laser-polarized Xe gas 
using NMR techniques {jj. 

Persistence has also remained a popular subject among 
applied mathematicians for many decades [^J. They are 
most interested in the probability of 'no zero crossing' of 
a Gaussian stationary process (GSP) between times T\ 
and T2 ||. It is well known that this probability usually 
decays as ~ exp(-0T) for large T = \T 2 - T t \ where 6 
is nontrivial ||||]. The persistence of some of the non- 
stationary processes mentioned in the previous paragraph 
such as the diffusion processes, can be mapped to that 
of a corresponding GSP [jlO|. This makes the two sets 
of problems related to each other and the power law ex- 
ponent in the former problem becomes the inverse decay 
rate in the latter. Even though 9 is, in general, hard 
to compute analytically, it is very easy to evaluate nu- 
merically in most cases. Given this fact, and the com- 
bined interest of both statistical physicists and applied 
mathematicians, much recent effort has been devoted to 
computing 9 numerically to extremely high precision. 

This raises a natural question: How accurately can 
one measure 91 Is there a natural limitation and if so, 
can it be overcome? This issue arises from the following 
simple observation. All the stochastic processes men- 
tioned above occur in continuous time. However, when 



one performs numerical simulations or experiments on 
persistence, one has to discretize time in some way and 
sample the data only at these discrete time points to 
check if the process has retained its sign. Due to this dis- 
cretization, some information is lost. For example, the 
process may have crossed and recrossed zero (or a spin 
may have flipped sign many times) between two consec- 
utive discrete time points. These crossings (or sign flips) 
go undetected due to the discrete sampling of the data. 
The question is how serious is this loss of information. Is 
it possible to estimate quantitatively the error involved 
due to the discretization? 

The purpose of this Letter is twofold: (i) to point out 
that there is indeed a very general and nontrivial effect, 
due to the discretization of time, on the measured per- 
sistence of any continuous stochastic process, and (ii) to 
provide a quantitative estimate of its magnitude in a sim- 
ple Markov model. The effect turns out to be nontrivial 
even for this simple toy model. We also develop two 
new analytical approaches, perturbative and variational, 
which provide results to extremely high precision. Wc 
emphasize that, even though we restrict ourselves here 
to a simple model by way of an example, this effect is 
very general and should be observable in simulations or 
experiments on more realistic systems. 

To formulate a precise quantitative question, let us 
consider a stationary stochastic process in continuous 
time T which is sampled at times T\, T2 , ■ ■ ■ , T n = T 
separated by a uniform window size, T — Tj_i = AT 
such that T = nAT. The continuous persistence P(T) 
is then approximated as P(T) w P n where P n is the 
probability that the process X(T) is positive at all the n 
discrete points. Note that, for finite AT, P n is different 
from P(T) since the process can cross zero more than 
once between two successive discrete times. One expects 
that the approximation P(T) « P n will improve as the 
window size, AT, decreases, and in the limit AT — > 0, 
n — > 00 keeping T = nAT fixed, P n — > P(T). By con- 
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trast, if the window size AT ^> r where r is the cor- 
relation time of the process, the stochastic variables at 
different discrete points become completely uncorrelated 
and we expect P n — > 2~ n , since the probability that at 
each point the process is positive is just 1/2. We then ask: 
How does the discrete persistence P n interpolate between 
these two limits as AT varies continuously from to oo? 
We show that for a GSP, in general, P n ~ [p{AT)] n for 
large n, where the function p(AT) is nontrivial with the 
limiting behavior 



p(AT) 



1 - 6 AT, 
1/2, 



AT 
AT 







(1) 



where 9 is the usual persistence exponent. As AT — > 0, 
one recovers the continuous persistence, P n — > (1 — 
6AT) n ~ exp(-0T) where T = nAT. The general goal 
would be to compute this function p(AT), the knowledge 
of which will provide an estimate of the difference, due 
to the finite window size AT, between the measured per- 
sistence P n and the P(T) of the underlying continuous 
process. 

The nonstationary processes discussed in the first para- 
graph are related to the equivalent stationary ones via 
T = hit Jl(|. A uniform spacing, AT, between mea- 
surements in the latter systems, therefore, corresponds 
to measurements uniformly spaced in log time in the for- 
mer. Such a measurement regime has indeed been used 
in a recent experimental study of diffusive persistence , 
with a spacing in log-time equivalent to AT w 0.24. The 
present paper is the first step in understanding how such 
discretization affects the measured result. To compare di- 
rectly with the experiment, we need to compute the func- 
tion p(AT) for the diffusion equation which is hard due 
to the non-Markovian nature of the process. However, 
to understand the general nature of this function p(AT), 
it would be useful to find a toy model where it can be 
computed explicitly. We consider below a simple Gaus- 
sian Markov process for which progress can be made in 
that direction. The physical process we study is the one- 
dimensional Ornstein-Uhlenbeck motion of a noisy, over- 
damped particle in a potential V(X) = pX 2 /2, where 
the position X of the particle evolves via the Langevin 
equation, 



dX 
~dT 



= -pX + V (T). 



(2) 



The white noise tj(T) has zero mean and a correlator 
( V (T)r,(T'))=2DS(T-T'). 

For this process, we first evaluate the continuous per- 
sistence and then compute the function p(AT). For the 
continuous persistence, a backward Fokker-Planck (BFP) 
approach is useful. Let Q(X,T) denote the probability 
that, starting at X at T = 0, the particle has not crossed 
the origin, X = 0, up to time T. We expect different be- 
havior depending on whether p > (stable potential) or 



p < (unstable potential). In the former case, the par- 
ticle will eventually cross the origin and hence Q(X,T) 
will decay exponentially with time. In the latter case, 
however, the particle has a finite probability to escape to 
infinity, and hence persistence should decay to a nonzero 
number. The latter case is also related to the problems 
of escape from metastable states studied before [[ll| . 
The probability Q(X,T) satisfies the BFP equation, 



dQ 
dT 



D 



if 



dX< 



pX 



dQ 
dX' 



(3) 



with boundary conditions Q(0,T) = and Q(oo, T) = 1 
for all T, and initial condition Q(X, 0) = 1 for all X > 0. 
The solution is 



Q(X, T) = Erf 



-iiT 



-.X 



(4) 



where D' = D/p and Erf [a;] is the error function. For 
p > 0, Q(X,T) becomes separable in X and T for large 
T, Q{X,T) ~ e-^ T X, and decays exponentially with T 
for fixed X. This gives the persistence exponent 9 = p. 
For p < 0, however, Q(X, T) approaches the steady state 
solution Q(X) = Erf {X/y/2\D'\) as T ^ oo. We also 
note from Eq. (|^) that the critical case p = corre- 
sponds to ordinary Brownian motion, and taking the 
limit p -> in Eq. one recovers the known result, 
Q(X,T) = Erf[X/V4DT], which decays as a power law, 
q\x,T) ~ X/VT, for large T. 

For later purposes, we will also need the Green's func- 
tion G(X2, Ti\X\, Tx), the probability that the particle 
starting at X — X\ at T = T\ will reach X-i at Ti, with 
T2 > Ti . This propagator can be easily computed exactly 
from Eq. (||) and we get, 



G(X 2 ,T 2 \X 1 ,T 1 ) 



(X 2 - 



^/2t:D'{1 



. g 2D'(l-a 2 ) 



(5) 



where a = g-M^-Ti) Note that for /i > 0, < a < 1, 
while for p < 0, a > 1 (and D' = D/p < 0). 

We now turn to the discrete persistence P n of the con- 
tinuous process in Eq. (||). Let Q n (X) be the probability 
that starting at X at T = 0, the process is positive at all 
the discrete points T x = AT, T 2 = 2 AT, . . ., T„ = nAT 
separated by the uniform window size AT. Then the 
discrete persistence is P n — L Q n (X)Po(X)dX , where 
Po(X) is the distribution of the initial position of the 
particle and can be arbitrary. Using the Markov prop- 
erty of the process in Eq. (^), it is easy to write down a 
recurrence relation for Q n (X), 

POO 

Q n+1 (X) = / G(Y,AT\X,Q)Q n (Y)dY, (6) 
Jo 

where G is the propagator as in Eq. (|5|) with a = e _/iAT 
and Qq(X) = 1 for all A > 0. This recurrence is the 
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discrete analogue of the continuous BFP equation (0). 
Indeed, it can be checked that Eq. (||) reduces to Eq. (||) 
in the limit AT — > 0. To simplify the algebra, we con- 
sider the rescaled variable, x = X/ y/D'(l — a 2 ), in terms 
of which the recursion reads 



1 



exp[-(y - ax) 2 /2]Q n (y)dy, (7) 



where we have used the explicit expression for G from 
Eq. (§). 

Let us first consider the case p > 0, i.e., < a — 
e -^iAT where, guided by the continuous case, we 

expect Q n (x) — > p n q(x) as n — > oo at any fixed x. Sub- 
stituting this asymptotic form into Eq. (Q), we get an 
integral-eigenvalue equation for q(x), 



pq{x) 



exp[-(y-ax) 2 /2}q{y)dy, (8) 



with eigenvalue p(a) that evidently depends continuously 
on a. Although Eq. (|J) admits many eigenvalues, we are 
interested only in the largest eigenvalue since it domi- 
nates the asymptotic behavior of Q n (x) for large n. We 
also note that Eq. (ph determines the eigenfunction q(x) 
only up to an overall multiplicative constant. Let us first 
consider the limit a — > or equivalently AT — > oo. In 
this case, Eq. (§|) can be solved exactly to give p = 1/2 
and q(x) — const, thus recovering the correct limit- 
ing behavior, Q n (x) — > const 2~™ (and the initial con- 
dition, Qq(x) = 1 for x > 0, fixes the constant at 
unity). For small a, by expanding Eq. (|8|) in a Taylor 
series, it is easy to compute p{a) perturbatively, giving 
p = h + —a + 0(a 2 ). The goal now is to evaluate p(a) for 
arbitrary a. To this end we develop below two analytical 
approaches and compare them with the direct numerical 
integration of Eq. (0). 

Perturbative approach: We expand the factor 
exp(axy), from the exponential in Eq. (|§|), as a power 
series and integrate term by term, to get 



pq(x) = 



cxp(— a 2 x 2 /2) ^ b„ 



E 7= ^ 



(9) 



n/2 fx 

b n = —l dyy n exp(~y 2 /2)q(y). (10) 



n\ Jo 



Substituting (|jj) into ( ]To| ) leads to the matrix eigenvalue 
equation 



pbj\ — ^ ^ A nrn ^T) 



(11) 



m=0 



-4, 



y/4ir(l + a 2 ) \1 + 



2a 



(n+m)/2 r i n+m+1 



A rh 1 - ^ 

ynlmi 



This approach converts an integral eigenvalue equation 
into a matrix eigenvalue equation, with matrix elements 



that decrease exponentially as n and m increase. Com- 
puting the largest eigenvalue of the N x N submatrix 
(n, m = 0, 1, . . . , N — 1) gives a rapidly converging series 
of estimates for p as N increases. For a given N, the re- 
sult is exact to order e" -1 , where e = 2a/(l + a 2 ). In this 
way one can easily obtain results for p(a) correct to one 
part in 10 12 . Convergence becomes progressively slower 
as a — > 1, which is expected since e — ► 1 in this limit. 
For a — * 1, however, we have the analytical result p — * a 
[such that p" — > exp(-npAT) — exp(—(iT)], since we 
must recover the continuum result in this limit. 

Variational approach: It is possible to derive a use- 
ful variational inequality for p. First we note that the 
integral operator in Eq. (^), asymmetric in x and y, 
can be made self-adjoint via the substitution, q(x) = 



g{x) exp 



(i-g ) -21 



which gives, 



K{x,y)g(y)dy, 



where i^(x, y) = if (y, x) = exp[— -^-^-i(a; 2 + y 2 ) + axy]. 
Let /(a:) be any normalizable function, f 2 (x)dx = 1. 
Using elementary properties of linear vector spaces and 
the self-adjoint property of the integral operator, it be- 
comes evident from Eq. (|l|) that the largest eigenvalue 
p satisfies the inequality, 



(13) 



> 



1 



OO f oo 



f{x)K{x 1 y)f{y)dxdy. 



(14) 



One can then use any trial function f(x) containing one 
or more variational parameters and then maximize the 
right hand side of Eq. (JlJ) with respect to these parame- 
ters to derive a rigorous lower bound for p(a) for arbitrary 
< a < 1. 

The limiting forms of the true eigenfunction g(x) in 
Eq. (13) for a — » and a — * 1 can be easily worked 
out, and suggest a trial function of the form f(x) — 
A(b + x) cxp(— \x 2 /2). The amplitude A is fixed by the 
normalization condition, f 2 {x)dx = 1, while b and A 
are the two variational parameters. The right-hand side 
of the inequality in Eq. (|l3) can then be evaluated in 
closed form and the optimization with respect to b and A 
performed. The resulting variational estimate turns out 
to be very accurate for all < a < 1, when compared 
to numerical results, and agrees with the perturbative 
results to at least 4 or 5 decimal places. 

Numerical Integration: It is not difficult to integrate 
Eq. ^ directly. However, since Q n (x) — > 1 as x — > oo, 
numerically it is convenient to first make the transfor- 
mation Q n (x) = G n (x) exp[(l — a 2 )x 2 /4] in Eq. (^) and 
then study the resulting equation for G n (x) by numerical 
iteration, with an arbitrary initial condition. For large n, 
G n (x) converges to p n g{x) where g{x) is the solution of 
Eq. (|l3|). The eigenvalue p is determined from the slope 
of the log- linear plot of A n — J °° G n {x)dx ~ p n versus n. 
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In Table 1, we compare the numerical, variational, and 
perturbative estimates of p. The differences are small in 
all cases, and the variational bound is satisfied. 

The eigenfunction q(x) of Eq. (||) can also be calculated 
by using the series ([)]), with the coefficients {b n } obtained 
from the corresponding eigenvector of the matrix A, Eq. 
(|l2|). It is shown, for a = 0.5, as the lower curve in 
Fig. 1. The asymptotic large- a: behavior (dashed curve) 
can be obtained analytically by noting that for large x 
we can set the lower limit in (11) to minus infinity with 
negligible error. The resulting equation can be solved 

exp(X 2 /4)A,(A), 



12| 



with solution q(x) 



exactly 

where X = (1 — a 2 ) 1 / 2 ^, D V {X) is the parabolic cylinder 
function, and v = lnp/ In a. The asymptotic behavior 
is q{x) ~ x v . The variational trial function, however, 
misses this asymptotic behavior (see Fig. 1) even though 
the variational eigenvalue is very accurate. 



a 


Pnum 


Pvar 


Ppert 


1.0 


l 


l 


l 


0.8 


0.8524547 


0.852440 


0.852454696506 


0.6 


0.7405959 


0.740589 


0.740595939159 


0.4 


0.6477666 


0.647765 


0.647766585747 


0.2 


0.5684903 


0.568490 


0.568490321623 


0.0 


1/2 


1/2 


1/2 


-0.2 


0.4408132 


0.440813 


0.440813209205 


-0.4 


0.3900580 


0.390004 


0.390057988652 


-0.6 


0.3469679 


0.346814 


0.346967773049 


-0.8 


0.3106439 


0.310444 


0.310643770245 


-1.0 


0.2800859 


0.279890 


0.280085758710 



Tablcl. Estimates of the eigenvalue p(a) for —1 < a < 1, 
from numerical, variational and perturbative methods. 
The latter is the most precise, being accurate to the num- 
ber of figures quoted. 

Although Eq. (0) was derived for a > 0, one can also 
study this equation or, equivalently, Eq. (||) and Eq. (13), 
for negative a. Is there a physical meaning for negative 
a? Let R n (x) denote discrete 'alternating' persistence, 
being the probability that, starting at x > (x is re- 
lated to X as before) at T = 0, the particle's position 
changes sign at alternate discrete points up to the n-th 
step. Then R n (x) evolves via the recurrence equation, 



1 f° 

R n+ i{x) = —=l exp[-(y 

V27T J-oo 



ax) 2 /2}R n (y)dy. (15) 



Changing y — > — y inside the integral, and using R n (y) = 
Rn(—y) (since the process has zero mean), we find Eq. 
( |l5| ) reduces to Eq. (m) with a replaced by —a. Thus, 
R n (x,a) = Q n {x, —a) and hence the largest eigenvalue 
p{a) for negative a governs the asymptotic decay of 'al- 
ternating' discrete persistence. We also note that while, 
for a > 0, Q n {x) ~ p n (a)q(x) for large n only for a < 1 
(for a > 1, Q n (x) approaches a steady state - see later), 
for negative a, Q n (x) ~ p n (a)q(x) for all a < 0. Fur- 



p(l/a) = \a\p(a), which can be used to obtain p (and the 
corresponding eigenfunction) for a < — 1 from the results 
for — 1 < a < 0. In particular, p — > 1/2 for a — > implies 
p — > l/2|a| for a — > — oo. 

Finally, we turn to the unstable potential, p < 0, i.e. 
a = e _AlAT > 1. As in the continuous case, we expect 
that the solution of Eq. (0) for a > 1 will reach a steady 
state for large n, Q n (x) — > q(x), where q(x) will satisfy 
Eq. (||), but with p = 1. Evidently q(x) will depend on 
a, and in the limit a — > 1 + (i.e. AT — * 0) it reduces to 
the continuous result obtained from Eq. (Q). For general 
a > 1, it is again possible to obtain accurate variational 
and very accurate perturbative estimates for q(x). We 
omit the details here since they are somewhat similar 
to the a < 1 case. In Fig. 1, we plot the perturbative 
q(x) for a = 2 (upper curve). The variational result, 
and the numerical result obtained from direct iteration 
of Eq. (|7]), are both indistinguishable from the plotted 
curve. Note that the case a < — 1, discussed in the previ- 
ous paragraph, corresponds to alternating persistence in 
an unstable potential, which does not approach a steady 
state. 




x, 5x 

FIG. 1. The eigenfunctions q(x) for a — 0.5 (lower curve) 
and a = 2.0 (upper curve, abscissa = 5a;). Solid lines - pertur- 
bative results; long-dashed - variational; dashed - asymptotic 
result q(x) ~ x v , with v — hip/ In a ~ 0.530661 for a = 0.5. 

In summary, we have shown that the discrete per- 
sistence due to the finite size of the time windows dif- 
fers considerably from the continuous persistence usually 
studied and we have computed explicitly this nontriv- 
ial effect analytically for a simple Markov model. The 
work extending some of the techniques developed here to 
more realistic non-Markov processes is in progress. We 
conclude by noting the recent examples of discrete time 
persistence in dynamical systems MM . 



thcrmorc, from Eq. (12), one has the symmetry relation 
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